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We use a quasi 3-D thermal evolution model for a spherical comet nucleus, 



. which takes into account the diurnal and latitudinal variation of the solar flux, 

^ . but neglects lateral heat conduction. We model the thermal evolution and activ- 

O , ity of Comet 9P/Tempel 1, in anticipation of the Deep Impact mission encounter 

. with the comet. We also investigate the possible outcome of a projectile impact, 

assuming that all the energy is absorbed as thermal energy. An interesting re- 
sult of this investigation, is that the estimated amount of dust ejected due to 



P ! the impact is equivalent to 2-2.6 days of activity, during "quiet" conditions, at 



c/3 ! perihelion 



We show that production rates of volatiles that are released in the interior of 
the nucleus depend strongly on the porous structure, in particular on the surface 
to volume ratio of the pores. We develop a more accurate model for calculating 
^ ! this parameter, based on a distribution of pore sizes, rather than a single, average 

pore size. 
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1. Introduction 

Comet 9P/Tcmpel 1 is the target of NASA's Deep Impact mission (Belton and A'Hearn 
1999, Meech et al. 2000, A'Hearn et al. 2005). The spacecraft, successfully launched on 
January 12, 2005, and due to encounter the comet on UT 5:52 (± 3 min) July 4, 2005, will 
collect images of the nucleus, and release a smaller projectile spacecraft. The projectile will 
maneuver its path for a collision with the nucleus of the target comet. 

Comet 9P/Tempel 1 was first discovered in 1867, and recovered in 1873 and 1879. In 
1881, the comet passed close to Jupiter, its orbit changed and it was lost. Later investigations 
(Marsden 1963) revealed that it had undergone two more close encounters with Jupiter, in 
1941 and 1957. The comet reappeared in 1967 and since that time it has been observed at 
every apparition. Its present orbital period is 5.5 years and the perihelion distance is 1.506 
AU. 

In anticipation of the mission, comet 9P/Tempel 1 has been intensively observed by 
both professional and amateur astronomers (Lamy et al. 2001, Meech 2002) and data about 
its size, rotation, production rates, etc., has accumulated (Fernandez et al. 2003, Belton et 
al. 2005). On the theoretical side, models have been developed to simulate the impact and 
its consequences (Schultz and Anderson 2005, Wiinnemann et al. 2005). It is impossible, 
however, to simulate the impact accurately, since very little is known about the properties 
of cometary material. Clearly, the impact energy — originally kinetic — will divide between 
mechanical and thermal energies, but the proportion will be strongly determined by the 
nature of cometary material: its strength, porosity, thermal conductivity, and so forth. 

Since the dynamical side of the impact has already been investigated in detail (Nolan 
et al. 1996, Housen 2002, Schultz et al. 2002), we focus in this paper on the thermal aspect. 
We first choose an initial model that matches as well as possible the observations of comet 
9P/Tempel 1. This is achieved by running a number of models, for different assumptions 
and parameter combinations, through several orbital revolutions. The composition of these 
models will be assumed to include water ice, dust and 5 additional volatile species. A 
"thermal" impact on this model will then be simulated. 

The evolution — both long-term, prior and following the impact, and short-term during 
the event itself — will be calculated by means of a quasi-3D code (Prialnik et al. 2005). 
This code takes into account diurnal and latitudinal variations, but neglects lateral heat 
conduction, and is an expanded version of the code used by Cohen et al. (2003). 

In Section 2 we briefiy describe the numerical code. Section 3 describes improvements 
implemented to the code that take into account a pore size distribution, rather than an 
average (fixed) pore size. Next, in Section 4 we describe ground-based observations that have 
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been recently carried out in anticipation of the Deep Impact mission in order to determine 
the comet's activity levels as a function of heliocentric distance. We include a subset of the 
observations in this paper for the purpose of modeling the dust coma so that we can obtain 
dust fluxes and the heliocentric distances between which the comet is active as constraints 
and tests for the thermal models. In the following Section 5 we show the results of numerical 
calculations for the thermal evolution of comet 9P/Tempel 1 along several orbital revolutions. 
The simulation of the impact and its consequences are described and discussed in Section 6. 
Our main conclusions may be found in Section 7. 



2. Numerical model and peirameters 

Our model assumes a porous, spherical and initially homogeneous nucleus composed of 
amorphous and crystalline water ice, dust, and 5 other volatile species: CO, CO2, HCN, 
NH3, and C2H2. The equations of energy and mass conservation for this system are briefly 
summarized below (using standard notation, e.g., Prialnik et al. 2005). 

We use a quasi-3D approach, where diurnal and latitudinal temperature variations are 
calculated as they result from uneven surface heating by solar radiation onto a spinning 
comet. Lateral heat conduction is neglected, and so heat is assumed to flow radially (per- 
pendicular to the surface). Thus, different points on the surface do not interact. This simple 
approach is justified by the low thermal conductivity of cometary ice mixtures and by the 
thinness of the skin-depth (see Cohen et al. 2003). 

The bulk mass density p and porosity ^ are given, respectively, by 

P = Pa + Pc + Pv + ^{Ps,a + Pg,a) + Pd, (1) 
a 

* = 1 - (Pa + Pc)/Qice - ^ Ps,al Qa " Pd/ Qd, (2) 

a 

where q is the specific density of a species a. The meaning of indices is: a - amorphous water 
ice, c - crystalline water ice, ice - when the phases are indistinguishable, v - water vapor, s 
and g - solid and gaseous phases of a volatile species, respectively, and d - dust. The partial 
pressure, assuming an ideal gas, is 

p _ T^gPg,aT 

Because amorphous ice has a tendency to convert to the crystalline form spontaneously, we 
have 

dpa 



dt 



= -Ap„, (4) 
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where A is the rate of crystallization. If / is the total fraction of occluded gas, fa — f, 
the equations of mass conservation are 



^ = (1 - /) Va - Qv, (5) 



^ + V-J, = g,, (6) 
for H2O in both phases, and similarly, 

^+y-J.^faXPa + qa: ^ ^ -Qa. (7) 

for the other volatile species. Taking into account energy conservation throughout the nu- 
cleus, and combining with the mass conservation equations, we obtain the heat diffusion 
equation in the form 



a \ a / a 



(8) 



where Tiac is the heat released upon crystallization, and Ha is the heat of sublimation. 
The above set of time-dependent equations is subject to constitutive relations: u{T), A(T'), 
qa{T, ^, Tp), Jq(T, ^, Tp), K{T, Tp), where Vp denotes pore radius. These relations require 
additional assumptions for modeling the structure of the nucleus. 

Measurements by Schmitt et al. (1989) have shown that the rate of crystallization is 
given by 

A(r) = 1.05 X lO^V^^^^/T s-^ (9) 

The rate of sublimation — mass per unit volume of cometary material per unit time — is 
given by 



(-Pvap,a(-^) -fa) 



2nRgT 



(10) 



where the term in square brackets represents the sublimation rate per unit surface area, and 
S represents the surface to volume ratio, which is a function of the given porosity and pore 
radius. The saturated vapor pressure, Pyap, is given by the Clausius - Clapeyron equation: 

Pv,p = Poexp-^/^, (11) 

where Pq = 3.56 x 10^^ ^ m-^, and B = 6141.667 K (Fanale and Salvail 1984). 

Dust is assumed to be, in part, dragged along with the gas flowing through pores and, 
in part, lifted off the nucleus surface by the sublimating vapor. For the former, the dust 
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velocity is assumed to be equal to the gas velocity (see Podolak and Prialnik 1996). An 
efficiency factor is calculated, to take account of a dust size distribution that allows only 
grains up to a critical size to be dragged or lifted off. The rest may accumulate to form a 
dust mantle. The efficiency factor may be adjusted so as to allow or prevent the formation 
of a sealing mantle (see also Prialnik et al. 2005 and references therein). 

The equations of mass conservation and energy transport arc second-order in space and 
hence each require two boundary conditions. One boundary condition for eq.(8) is vanishing 
heat flux at the center. The second one is obtained from the requirement of energy balance 
at the surface: 

F(R) = .aT(R, tr + TP^iT) ^f^H - (1 - A)^^,cos (12) 

In the simple case considered here, of a rotational axis that is perpendicular to the orbital 
plane, the solar zenith angle is given in terms of latitude ^, and hour angle </? = (27r/Prot)^) 

cos 2; = cos6' cos 99. (13) 

The factor T < \ represents the fractional area of exposed ice, since the surface material is 
a mixture of ice and dust, and can be written as 

jr^ A ^ ftccAi\ ' ^^^^ 

V Pice Qd) 

(Crifo and Rodionov 1997). 

Similarly to the heat flux, the mass (gas) fluxes vanish at the center. At the surface, i?, 
the gas pressures are those exerted by the coma; in the lowest approximation they may be as- 
sumed to vanish: Pa{R, t) = 0. It should be mentioned that the simple (and commonly used) 
outer boundary conditions for both energy and gas fluxes have been recently examined in 
more detail by Davidsson and Skorov (2002, 2004). Back-scattering of molecules leaving the 
nucleus surface and penetration of solar radiation into a thin sub-surface layer of the nucleus 
have been shown to affect the surface and sub-surface temperatures. These temperatures, 
however, are equally affected by the other approximations of the model (such as sphericity 
or homogeneity of surface structure) , not to mention the uncertainty in the thermal conduc- 
tion coefficients. Nevertheless, production rates are far less affected by these factors. For 
example, if back scattering reduces the net amount of sublimation at a given temperature, 
the surface temperature increases, but with it the sublimation rate increases as well. In fact, 
at low heliocentric distances, the amount of sublimation may be quite accurately calculated 
simply by (1 — A)Lq cos z / {AikP^H,) , independently of the surface temperature. 
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Table 1: Orbital elements and nucleus characteristics 



Parameter 


Symbol 


Value 


Units 


Semi-major axis 


a 


3.12 


AU 


Eccentricity 


e 


0.517491 




Effective radius 


R 


3.3 


km 


Geometric albedo 


A 


0.04 




Nucleus spin period 


Plot 


41.85 


hr 


Dust mass fraction 




0.5 




Ice mass fraction 


^ice 


0.5 




Porosity 


^ 


0.5 




Dust density 




3250 


kg m~^ 


H2O ice density 


^ice 


917 


kg m~^ 


Dust heat capacity 


Cd 


1.3 X 10^ 


J kg-i K-i 


Ice heat capacity 


Cice 


7.49T + 90 


J kg-i K-i 


Dust conductivity 




10 


J m-i s-i K-i 


C-Ice conductivity 




5.67 X IOVT 


J m-i s"^ 


A-Ice diffusivity 


■^a./ ( Cice ^ice ) 


3 X 10"^ 


2 —1 
m s 



The numerical code solves the heat diffusion and mass conservation equations using a 
fully implicit difference scheme (see Prialnik 1992). In order to obtain a better resolution 
near the surface, where the temperature gradient is steepest, the layers are progressively 
thinner towards the surface. With the spin vector perpendicular to the orbital plane, results 
for 6 hour angles are recorded: (/? = 0/360°, 60°, 120°, 180°, 240°, 300°, a different one at 
each time-step. Separate calculations are carried out for latitudes corresponding to cos 6 = 
1, 0.75, 0.5, 0.25, that is, between the equator {6 = 0°) and a near-pole angle {9 = 75.5°). 
The physical parameter values used in the model calculations are given in Table 1. 

3. Properties of a porous nucleus structure 

A cold, icy porous medium allows for sublimation from the pore walls and condensation 
onto them, as well as for flow of — usually dilute — gas through the pores. The three 
fundamental properties that affect these processes are the porosity the surface to volume 
ratio (SVR) S and the permeability 0. The first is defined as the total volume of pores per 
given bulk volume, and the second, as the total interstitial surface area of the pores per given 
bulk volume. The permeability is, up to a numerical constant, the proportionality coefficient 
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between the mass flux and the gradient of P/ \/T that drives the flow of mass, where P and 
T are the gas pressure and temperature, respectively. 

In order to determine these three properties, a model of the porous structure is required. 
A commonly used one in modeling cometary interiors is that of a bundle of unconnected 
cyhndrical tortuous capillary tubes (Mekler et al. 1990), where the tortuosity ^, defined as 
the ratio of capillary length to linear distance, is taken as a free parameter. Another free 

parameter is the pore (capillary) radius r, for which some reasonable average pore size is 
assumed. These assumptions result in very simple expressions for S and (j) in terms of ^, ^ 
and r: 

2\1/ 

S=- 0=^. (15) 

These expressions imply, however, that all capillaries are of the same radius. If r represents 
the mean of a distribution of pore sizes, the simple formulae above cease to be correct, a fact 
that has been largely overlooked in previous calculations. There is, however, a price to pay 
for a more realistic approach to pore sizes, and it involves a larger number of parameters for 
defining the medium. Nevertheless, we shall show that the deviations from eqs.(15) may be 
quite considerable and thus worth the price. 

Starting from the same model of cylindrical tortuous capillaries, we assume the radii 
to vary according to some size distribution, and consider a volume of unit thickness. Let 
N[r)dr be the number of capillaries with radii between r and r + dr crossing a unit area. 
Keeping in mind that for a capillary tube, (p oc r^/^ (Gombosi 1994), the three fundamental 
properties of a porous medium are, in fact, moments of the distribution function: 

S ^ i j 27rrN{r)dr (16) 

^ = ^ /" 7rr^N{r)dr (17) 



= I y r^N{r)dr (18) 



(cf. Prialnik et al. 2005). It follows that 

e '!r''N{r)dr ~ C 



* r r^N(r)dr * , , 



S = 2^ i )\ = 2^ - , (20) 



while 

Q — 9\Tj ^ ^ 

J r'^N{r)dr \r ^ 

where f in eq.(19) is the mean pore radius weighted by the volume fraction occupied by 
capillaries of radii in the range (r, r + dr). We note that eqs.(19) and (20) are similar to 
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(15), except that S depends on the harmonic mean of (^), rather than on the inverse of the 
mean. If we adopt f as the free parameter, then the expression for S in eq.(15) should be 
corrected by a factor 

''\r) {Jrm{r)drY 



A reasonable distribution function for pore size is a power law, since grains ejected 
from comet nuclei have a power law size distribution (sec below), and comets are believed 
to be formed by uncompacted aggregations of grains, in which case the grain and pore size 
distributions should be similar. Thus, we assume 

N{r) oc r"'* for r^in < r < r^ax, (22) 

which depends on 4 parameters: the normalization factor, the exponent a and the range 
of the distribution defined by rmin and r^iax- Instead, besides a, we may choose to assume 
values for the porosity ^ (serving as a normalization factor), the average pore radius f (in 
order to be able to compare results with the former, simpler approach), and the ratio 

X = Tmin/rmax < 1- (23) 

The resulting correction factor for the SVR, 

r (r.\ (3-«)^ (i-x^-")(i-x^-") 

is plotted in Fig.l for several values of X. Following the distribution of dust grains inferred 
from observations, typical parameter values are X = 10~^ — 10~^ and 3 < a < 4 (FuUe et 
al. 1997, Jockers 1999, Barker et al. 2002). 

We note, therefore, that for given porosity and average pore size, taking into account 
a pore size distribution can increase the SVR by up to a factor of 100. Moreover, large 
corrections correspond precisely to that narrow range of a values deduced from observations. 
The effect of a large SVR should manifest itself both externally, through the production 
rates of volatiles, and internally, through composition profiles. Both are bound to affect, 
for example, the interpretation of the upcoming Deep Impact observations. A comparison of 
model results obtained with and without the correction factor for the SVR will be shown in 
Section 5. 



4. New Ground-Based Observations 



In order to prepare for the Deep Impact encounter, the Deep Impact team has undertaken a 
large observing campaign to characterize the nucleus and the levels of activity (Meech et al. 




Fig. 1. — The correction factor of the SVR parameter, as a function of the power law expo- 
nent, for different ratios of minimal to maximal pore size. Note that the interval commonly 
used for pore size distribution (3 < a < 4) gives a factor > 10. 
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2005). Prom past apparitions it was known that this comet typically exhibits a sharp rise in 
brightness near 200 days pre-perihelion (Meech et al. 2005). Optical CCD images have been 
obtained regularly since 1997 to monitor the development of activity and the cessation of 
activity as the comet moved to aphelion in early 2002. This comet is well placed for observing 
every other apparition, and the early 2000 perihelion passage was only moderately good for 
dust dynamical modeling, but was our first chance to measure the dust parameters fot this 
comet. Nevertheless, we have a large data set accumulated for this purpose beginning in 
January and March 1999 pre-perihelion, and continuing from August 2000 through January 
2001. 

Finson-Probstein (1968) dust-dynamical models use the observed extent and morphol- 
ogy of the dust coma to determine the relative velocity distribution, size distribution and 
production rates of the dust leaving the nucleus as a function of heliocentric distance. These 
models evaluate the motion of a suite of particles after leaving the nucleus under the in- 
fluence of solar radiation pressure and gravity. The scattered light from the dust is added 
together and fit to the surface brightness of the observed coma. 

In order to use dust imaging to constrain the velocity and size distributions as well as 
the production rates, the images must be taken at the appropriate times: 

• Small particles move rapidly away from the nucleus (and out of the field of view) , so 
in order to constrain the size distribution from models, many images closely spaced in 
time need to be obtained. 

• Large grains, which move slowly along the orbit, need observations equally spaced over 
time periods of months for proper modeling. 

• Small particles tend to lie in the anti-solar direction, whereas large particles lag behind 
in orbit. 

• The predicted dust trajectories (syn-curves, comprising the locus of the synchrones 
and syndynes) should be spread out in the plane of the sky so that good constraints 
can be obtained for the model parameters. 

• Very circular orbits tend to have widely spaced syn-curves, but changes with time are 
harder to see. Comets on very elliptical orbits have problems inbound when the large 
particles fall along the sun-to-comet radius vector (and overlap the syn-curves for the 
small grains). 9P/Tempel 1 is a good comet for modeling in that its orbit has an 
intermediate eccentricity. 



For comet 9P/Tempel 1, the syn-curves were predicted to be moderately good for the 
fall of 2000 observations (post-perihelion), and excellent starting in February 2005 through 
the encounter. We will report on the preliminary modeling for the observations obtained 
during fall 2000. 

Kron-Cousins R-band images were obtained on 2000 Aug 21 and 2000 Sep 30, using 

the University of Hawaii 2.2m telescope on Mauna Kea and the the Tektronics 2048x2048 
CCD (read noise = 7e~, gain = 1.4e^/ADU, plate scale = 0.219'.' pixel^^). Flat fields were 
obtained on the twilight sky, and standard reduction procedures were able to flatten the 
data to with 0.4% across the CCD. The nights were photometric, and standard stars from 
Landolt (1992) were observed on each night. 

Composite images were made by using measurements of the centroids of a large number 
of field stars and the centroids were used to compute the offsets between the individual 
images (from telescope guiding errors). After the offsets were applied, the ephemeris rates 
for 9P/Tempel 1 were used, in combination with the image plate scale, to calculate shifts to 
simulate guiding at non-sidereal rates, and then the images were added. Selected field stars, 
close to the comet, were fit and removed from composite image. The resultant images are 
shown in Figure xx (attached to the paper). 

Fig.xx.- Composite images of comet 9P/Tempel 1 (left) from 2000 Aug. 21 (r=2.54 AU), 
composed of 26 x 300 sec R images and (right) from 2000 Sep. 30 (r=2.77 AU) composed 
of 42 X 300 sec R images. Images are 180x180", with N at the top and E to the left. 

Using a code developed by Farnham (1996), Finson-Probstein dust dynamical models 
were fit to the images, and preliminary results are shown in Table 2. Contour fits to the 
images are shown in Figure 2. 



Table 2: Finson-Probstein Model Fit 



Parameter 



Value 



Smallest particles 
Largest particles 
Emission start 
Maximum dust output 
Emission dechne 
Velocity {v{(3))^^^ 



a 



[2] 



3 jim. 
3 mm 

perihelion - 100 days 
perihelion - 60 days 
perihelion -\- 260 days 



u(r„ 
3.1 



i)/'y(rmax) ~ 40 



Notes: ^^^fi = Frad/ Fyrav = 5.74 x 10 ^Qrp/Qdi", is a proxy for grain size, where Q is the radi- 
ation pressure scattering efficiency. Fit value for the slope of the particle size distribution. 



Fig. 2. — Finson-Probstein dust dynamical model fits to the data on Aug. 21 (left) and Sep. 
30 (right), 2000 plotted on contours from the composite images. 

5. Results of evolutionary calculations 

Since little is known about any comet's structure and composition, our first task was to 
choose the structural and compositional parameters. The composition was chosen by running 
several models of various generic composition, but identical structure — as listed in Table 
3 — and comparing the resulting production rates to available observations (Cochran et 
al. 1992, Osip et al. 1992). This is shown in Fig.3. The model for which agreement was 
obtained for all observed species available was then chosen as the working model. The dust 
production rate of this model, shown in Fig.4 agrees very well with the results derived from 
ground-based observations described in Section 4. We note, in particular, the times of rise, 
maximum and decline in the rate of dust emission, compared to those of Table 2. The 
magnitude corresponding to the dust production rate at maximum is about 9; varying the 
average grain size and/or dust grain density may result in a change of ±0.5 mag (for a 
correlation between dust production and magnitude, and perihelion magnitudes, see Meech 
et al. 2005, 2002). 

Guided by the dust grain size distribution derived from observations, we adopted for 
the pore size distribution inside the nucleus a slightly steeper power law, with a = 3.5, an 
average pore size of 100 fim, and a ratio X = 10^^. According to the discussion in Section 3, 
^ ~ y/f\ninTmax, wMch rcsults in a size range between rmin =1 A'-ni and rmax =1 cm. The range 
was chosen to be wider than the dust grain size range (3 /xm-3 mm, see Table 2) in order 
to account for micro-pores at the low-end, and to enable the flow of a few mm-size grains 
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Table 3: Initial composition of 9P/Tempell models 







Ti'fi'n'npn cfat; 


Tpp mivfiiTP 


TpP mTYfllTP 


r VI Q cc 




T)nt;"t" Tn;^ri"t"lp 


IMo TTlPTI'l'lp 

1 >i ^ liidjiiijiC' 


T)nG;"t" TTli^TlI'lp 


fraction 


No mantle 










0.5 


0.5 
















0.45 


0.45 


Xd 


0.5 


0.5 


0.5 


0.5 


Xco 


0.05 


0.05 


0.025 


0.025 


Xcoi 


0.0125 


0.0125 


0.00625 


0.00625 


Xhcn 


0.0035 


0.0125 


0.00625 


0.00625 


XnHs 


0.0125 


0.0125 


0.00625 


0.00625 


XC2H2 


0.0012 


0.0125 


0.00625 


0.00625 



HjO production rate, on inbound leg (8th orbit taken for all models) 




+ Osipetal. 1992 

— Working model 

— Trap, gas, mantle 

— Ice mix., no mantle 
Ice mix., mantle 



HCN production rate, on inbound leg (8th orbit taken for all models) 













+ Osip et al. 1992 
o Cocliran et al. 1 992 












Working model 
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C.H production rate, on inbound leg (8th orbit taken for all models) 



1.9 



1.85 



1.75 



1.65 



1.6 



1.55 



1.5 



o Cocliran et al. 1992 

— Working model 

— Trap, gas, mantle 

— Ice mix., no mantle 
Ice mix., mantle 



Fig. 3. — Comparison of observation and model run results for 4 generic models, as described 
in the paper. Observed production rate of OH, as the product of H2O, CN, as the product 
of HCN, and C2, as the product of C2H2, were taken from the references. 
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Dust production rate (8th orbit) 

6.4 1 1 1 1 1 1 1 1 

6.2- 




Days (from perihelion) 

Fig. 4. — Dust production rate pre- and post-perihelion for the basehne model. 

requiring somewhat larger pores at the high-end. A steeper power law was chosen for pore 
sizes than indicated by grain observations, because dust grains are also lifted off the surface, 
where pore size does not constitute an impediment, and hence the relative proportion of 
larger grains coming from the surface should be larger than for those originating from the 
interior. 

Assuming the spin axis to be perpendicular to the orbital plane and starting with a 
uniform low temperature, we followed the evolution of the working model for several orbits. 
The thermal and compositional evolution is illustrated in Fig. 5, for the subsolar point on 
the equator and for a point near the pole. These represent extreme cases — maximum and 
minimum — with respect to absorption of solar radiation, and hence activity, of a spinning 
nucleus in the orbit of 9P/Tempel 1, among all possible inclinations of the spin axis. We 
recall that the model nucleus is spherical, whereas in reality the nucleus is found to be 
highly elongated, with semi-axes a ~ 7.2 km, 6 ~ c ~ 2.2 km (Belton et al. 2005). Thus the 
"collecting area" of the nucleus may vary, according to inclination, between ~ 15 — 50 km^, 
compared to 34 km^ of the model. This would place an error bar of up to a factor of 2 on 
the results regarding production rates, since most of the solar heat is spent in sublimation 
of volatiles. The shape itself should be of lesser importance, since the affected layer, as we 
shall see, is much smaller than the radius. 
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The outstanding feature emerging from this calculation is the complicated stratification 
pattern as a function of depth, where layers enriched in various volatiles alternate. Moreover, 
several layers enriched in the same volatile may appear at different depths. The effect is 
illustrated by the mass fraction of amorphous ice shown in the right panels of Fig. 5. We 
recall that the model's composition includes 5 volatile species trapped in amorphous water 
ice. These volatiles cover a wide range of sublimation temperatures (see Prialnik et al. 2005). 
As the surface of the nucleus is heated and the heat wave propagates inward, the amorphous 
ice crystallizes and the trapped gas is released. The gas pressure in the pores peaks at 
the crystallization front. As a result, gas flows in part outward and escapes, and in part, 
inward into colder regions. Eventually, each species reaches a sufficiently cold region for it to 
recondense. Since recondcnsation releases heat, it affects the composition of its surroundings 
and thus a complicated pattern results, of alternating ices mixed with the amorphous water 
ice. When another heat wave reaches these regions, on a subsequent perihelion passage, the 
heat is absorbed in sublimation of the recondcnscd volatiles rather than in crystallization of 
amorphous ice. This is how alternating layers of crystalline and amorphous ice arise, rather 
than a single boundary between a crystalline exterior and an amorphous interior. Finally, 
the steps that appear in the outer crystalline/amorphous ice boundary are due to erosion of 
the nucleus, which brings this boundary closer to the surface. 

The stratified layer extends from a depth of about 10 m below the suface and down 
to a few hundred meters. Since 10 m is roughly the orbital skin depth for 9P/Tempel 1, 
this structure should cause activity variations on the orbital time scale. This means that 
activity may differ from orbit to orbit and occasional spurious outbursts may arise, when 
the heat wave propagating down from the surface reaches a region enriched in ice of some 
volatile species. Such outbursts of gas should be accompanied by ejection of dust. Indeed, 
this variable behavior is exhibited in Fig.6, where we note an outburst of gas (CO and CO2) 
and dust following several "quiet" regular orbits. By contrast, water production, whose main 
source is sublimation from the surface, follows a much more regular pattern. 

The evolution of local noon temperature profiles is shown in the left panels of Fig.5. 
It indicates that a steady state is reached in each case after a few revolutions at about 
the same depth, both near the pole and at the equator. This depth corresponds to the 
orbital skin depth \2Ka^/'^ / {^/GM^pc)]^/'^ ^ 10 m. We note, however, that the maximum 
temperatures at the two locations differ by about 40 K; accordingly, the skin depth, which is 
roughly inversely proportional to temperature, is slightly larger near the pole. Without the 
correction factor to the SVR, internal temperatures are somewhat higher in the outer layers 
since less heat is absorbed in ice sublimation. 

We now wish to draw attention to the effect of the corrected SVR on the results. 
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Fig. 5. — Evolution of the comet nucleus working model through repeated revolutions. The 
temperature and amorphous ice abundance are shown as a function of time and depth: (a) 
and (b) — near-pole latitude, with the corrected SVR term; (c) and (d) — equator, with the 
corrected SVR term; (e) and (f ) — equator, with the original SVR term. The depth scale is 
logarithmic, from the surface to the interior. 
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Temperature profiles at the equator are compared in Fig. 5 in the two lower left panels, and 
the evolution of crystallization in the two lower right panels. The correction factor increases 
the internal pore surface and enhances sublimation from the pore walls. This hinders the 
heat wave from penetrating deeper and internal temperatures are thus lower. As a result, 
the depth of amorphous ice is shallower, as illustrated in Fig.5. At higher latitudes, due to 
the diminished insolation, it is still shallower. 

Production rates obtained for long-term evolution with and without the SVR correction 
are compared in Fig. 6. As expected, surface properties, such as temperature and H2O 
production rate, arc not affected. The production rates of gases released in the interior, 
either escape from crystallizing ice or by sublimation of rccondcnsed ice, are affected to a 
larger extent, but preserve the same pattern of behavior. This is not surprising, since the 
driving energy source is the same. We note that the complex stratified structure of the 
outer layer of the nucleus gives rise to occasional outbursts. Clearly, the orbital evolution 
of production rates differs considerably among different volatiles. This means that observed 
volatile abundances do not reflect nucleus abundances (cf. Huebner and Benkhoff 1999, 
Prialnik et al. 2005). 



6. Simulation of a deep impact 

The aim of the Deep Impact mission is the collision of the impactor spacecraft with the nu- 
cleus, in order to create a crater. This will enable the flyby spacecraft to perform observations 
and measurements of the impact itself, the ejected material and the interior composition, as 
revealed by the exposed crater. 

The thermal effect of the impact is simulated by assuming an additional energy influx 
for an infinitesimal period of time. The known input parameters for the impact simulation 
arc the total (kinetic) energy of the projectile -Etot) the estimated area A, that will be affected 
by the impact — which is needed in order to obtain an energy flux — and the heliocentric 
distance (prcpcrihelion) where the impact is expected to take place, d^. The total energy 
absorbed per unit area is Ea = Etot/A. We note that the relevant parameter is Ea, hence 
the calculation is not sensitive to the total impact energy or the area, separately 

Since numerical calculations do not allow singularities, we assume the energy deposition 
rate to be in the form of a narrow time-dependent Gaussian (rather than a delta function), 
centered at the time of the impact, with a width corresponding to the timescale for the 
formation of the crater, 

F(t) = Foe-'^l^*-*")/^!'. (25) 
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Fig. 6. — Comparison of original SVR (blue) and corrected SVR (red), for 10 orbits of the 
working model: results are shown for the surface temperature at the subsolar point as well 
as production rates for all volatiles and dust. 
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Table 4: Impact simulation properties 



Parameter 


Symbol 


Value 


Units 


Total energy 




1.9 X 10^° 


J 


Crater area 


A 


1000 


m^ 


Crater depth 


AL 


30 


m 


Heliocentric distance 


do 


2 


AU 


Revolutions before impact 


n 


7 




Timescale 


T 


180 


s 


Peak energy flux 


Fo 


10^ 


J m~^ s"^ 



Here r is a free numerical parameter that may be interpreted as the impact timescale. 
Normalizing the deposition rate, 

/oo 
F{t)dt = Ea (26) 
-oo 

we obtain Fq — Ea/t. Thus, the amplitude of the energy deposition rate is related to 
the physical parameters of the impact. The parameter values used are listed in Table 4. 
The resulting peak energy flux is about 300 times higher than the solar energy flux at that 
distance at the subsolar point. At a distance of 1.5 AU, the current impact distance, the 
contrast would be reduced by a factor of ~ 0.6, but this change should not be signiflcant. 
The time of impact to is calculated by 

to = V«V^M0(27rn-7r-'(?o + esin^?o) (27) 
i?o = cos-^[(l-V«)/e], (28) 

where n is the number of orbital revolutions calculated prior to the impact. 

At the onset of the impact — deflned to occur at a time tg ~ ^'T-, where v > 1, so that 
F{t ± vt) -C — we remove from the model nucleus a layer of thickness AL k, 30 m, to 
simulate the crater that is expected to form. Thus the thermal energy is deposited below 
the original nucleus surface, where the composition had been either preserved, or altered by 
earher evolution to a different extent. We should bear in mind that the thermal effect of the 
impact is maximized, since it is assumed that the entire kinetic energy is turned into heat 
and the heat is all deposited at a depth corresponding to the bottom of the impact crater. 

Two impact calculations were carried out for two different locations on the nucleus 
surface: one at the equator, and another near the pole. The purpose was to find out to 
what extent is the point of impact expected to affect the outcome. The comparison of the 
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two locations is presented in Fig. 7, and shows that, to within a fraction of a magnitude, the 
production rates of the major constituents (H2O, dust, CO and CO2) are indistinguishable 
at the moment of impact. However, since the flyby spacecraft is supposed to conduct the 
observation for a period of ~ 15 minutes after the impact, the emission of dust, CO and 
CO 2 may vary by ~ 2 orders of magnitude during this time. 

The low production rates before impact are an artifact of our simulation: when we 
remove an outer layer of several tens of meters, we expose cold material. Since the energy 
is supplied slowly at first, production rates drop. The total amount of dust ejected upon 
impact amounts to ~ 1.7 x 10^ tons at the ncar-polc location, and ~ 2.2 x 10^ tons near 
the equator. Compared with the dust production rate over the entire comet at perihelion, 
these amounts are equivalent to the total dust output of 2 days and 2.6 days, respectively, 
at perihelion. 

We note that the CO and CO2 ejection rates are not smooth: they show two spikes 
and a somewhat later peak (see Fig. 7). This is due to the layered structure of the nucleus 
described in Section 5. The fact that the behavior pattern of these volatiles is reflected in 
the rate of dust ejection indicates that a significant fraction of the dust originates in layers 
beneath the surface. Small dust grains are thus dragged by volatiles through the pores. This 
effect is less marked near the pole, where the subsurface layers have been processed to a far 
lesser extent during evolution prior to the impact. 

An interesting result of the long-term evolution, is that even though the thermal energy 
influx due to the impact is confined in both time (by the Gaussian function) and space (by 
the depth of penetration), the effect lingers and shows deviations in production rates - as 
compared with the undisturbed model - long after the event. However, this effect cannot be 
used to characterize a comet that has undergone an "impact", because the evolution pattern 
is similar to that of a comet model that was not disturbed. The difference is in the detailed 
variation of production rates with heliocentric distance, but these details depend to a similar 
extent on model parameters and assumptions. 

7. Conclusions 

We have followed the evolution of a rotating nucleus model in the orbit of comet 9P/Tempel 
1 for several orbits, until a quasi-steady-state was attained, and then we have simulated 
an impact, similar in energy to the projectile of the Deep Impact mission, focusing on the 
thermal outcome. A similar effect may be expected of a random collision with a large meteor. 
Provided that such a collision will create a deep crater and expose layers a few tens of meters 
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Impact at a near-pole latitude Impact at the equator 
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Fig. 7. — Comparison between the two extreme of impact latitudes: near-pole (75.5°) and 
equator (0°). The dominant volatiles, H2O, CO and CO2, and the dust component, are 
presented at high temporal resolution, around the peak of the impact. Note that the general 
shape of production rates is similar at both impact sites. 
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deep, an outburst of gas and dust ejection is expected to result. This is due to the fact that 
at such depths, ices of volatile species should be present, as well as amorphous water ice, 
which may crystallize by absorbing the impact heat and release more volatiles. Our main 
results regarding the outcome of the impact are as follows. 

• The increase in production rate of volatiles and dust is of several orders of magnitude, 
and thus readily observable. The total dust output, for example, is equivalent to the 
output of about two days at perihelion. We should bear in mind, however, that these 
results provide an upper limit, since it is assumed that the entire energy is spent in 
thermal effects. 

• An important conclusion of these calculations is that the place on the nucleus where 
the impact occurs is not as significant as one might expect. The total dust output, for 
example, differs by only a factor of ~ 1.3. 

• Variability in the comet's activity may be detected even on the short time scale of 
the impact, as the heat wave propagates into a layer of stratified composition, where 
amorphous water ice may crystallize and ices of various volatile species may sublimate. 

• Subsequent long-term thermal evolution is also affected to some extent, but not in a way 
that would be recognized by observations as an aftermath of an impact (or collision). 
The variation of volatile production rates with time (or heliocentric distance) will retain 
a general uneven pattern, but whether this pattern differs from that of an undisturbed 
comet may be difficult to establish. 

It should be borne in mind that these predictions are based on a simplified model of a real 
comet. 
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